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SUMMARY 


This paper presents an hp-adaptive discontinuous Galerkin method for linear hyperbolic conser- 
vation laws. A priori and a posteriori error estimates are derived in mesh-dependent norms which 
reflect the dependence of the approximate solution on the element size (h) and the degree (p) of 
the local polynomial approximation. The a posteriori error estimate, based on the element residual 
method, provides bounds on the actual global error in the approximate solution. The adaptive 
strategy is designed to deliver an approximate solution with the specified level of error in three 
steps. The a posteriori estimate is used to assess the accuracy of a given approximate solution and 
the a priori estimate is used to predict the mesh refinements and polynomial enrichment needed 
to deliver the desired solution. Numerical examples demonstrate the reliability of the a posteriori 
error estimates and the effectiveness of the hp-adaptive strategy. 


INTRODUCTION 

Adaptive methods, based on successive refinement of an existing mesh or a complete re-meshing 
of the computational domain, have become invaluable tools in computational fluid dynamics. The 
amount of refinement or the clustering of grid points is often determined by an element refinement 
indicator, 0 K , of the form 

e K = &“|D"u| (i) 

where h is a measure of the element size, a is an exponent, and D n u represents some higher-order 
derivative of a key variable. Since these indicators are based on interpolation or truncation error 
estimates, they are applicable to a large class of problems and are independent of the numerical 
method used to obtain an approximate solution. Although these indicators detect certain flow 
features, they may not relate to the actual error in the solution. Often, these indicators only 
provide a relative measure of the error and do not provide a criteria for stopping the adaptive 
process. 

This paper summarizes some of the work presented in [4] aimed at developing hp - adaptive 
strategies based on reliable error estimates for hyperbolic conservation laws. The work focuses on 
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the discontinuous Galerkin method applied to a model class of linear hyperbolic conservation laws 
for which it is possible to develop mathematically rigorous a priori and a posteriori error estimates. 

The notion of discontinuous Galerkin methods for hyperbolic problems originated in the classical 
work of Lesaint and Raviart [1] over two decades ago. Johnson and Pitkaranta [2] generalized the 
theory of discontinuous Galerkin methods by introducing mesh-dependent norms and were able to 
derive a priori error estimates in such norms for linear hyperbolic problems. Discontinuous Galerkin 
methods were extended to nonlinear hyperbolic conservation laws by Cockburn, Hou, and Shu [3] 
who developed a local projection strategy to provide nonlinear stability. 

The local nature of the discontinuous Galerkin method makes it ideally suited for adaptive 
strategies which combine local mesh refinement (h) with local enrichment of the polynomial ap- 
proximation (p) in an element to improve solution accuracy. The theory of discontinuous Galerkin 
methods was extended to hp-finite element approximations by Bey in [4] for a class of linear hy- 
perbolic conservation laws. In [4], very high accuracies and convergence rates were observed in 
applying discontinuous Galerkin methods to representative test problems. 

The a posteriori error estimates developed in [4] are based on the element residual method and 
provide bounds on the global error. Error estimates are combined with an /ip-adaptive strategy 
that predicts the mesh required to deliver a solution with the specified level of error. 

The theoretical developments of [4] are summarized in this paper. The discontinuous Galerkin 
method, the a priori error estimate which establishes the accuracy and convergence of the method, 
and the a posteriori error estimate used to assess the accuracy of the numerical solution are pre- 
sented. The reliability of the a posteriori error estimate is assessed by solving two examples prob- 
lems with known discontinuous solutions. The effectiveness of the adaptive strategy at delivering a 
solution with the specified level of error is" also demonstrated using the numerical examples. 


THE DISCONTINUOUS GALERKIN METHOD 


Consider the following hyperbolic conservation law 

f3 • Vm + au = / in OcU 2 (2) 

f3 ■ n u = (3 • n g on T_ (3) 

where (3 = (/3i,/ 3 2 ) T denotes a constant unit velocity vector, n denotes the unit normal vector 
pointing outward to the domain boundary <90, T_ = {x G dQ | j3 ■ n(x) < 0} denotes the inflow 
boundary, a = a(x ) is a bounded measurable function on 0 such that 0 < ao < o(x), / G L 2 ( fl), 
and g G L 2 (T_). While this is the simplest of hyperbolic conservation laws, solutions to (2) may 
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contain discontinuities along characteristic lines x(s) defined by = /3. Solutions to (2) belong to 
the space of functions V(Ct) = {«£ L 2 ( fl) | v@ € L 2 { 0)} where vp — (3 • Vv. 


The starting point for the discontinuous Galerkin methods is to develop an appropriate weak 
formulation of (2) defined on a partition of fl into elements, denoted by Vh- Here the elements 
K G Vh are general quadrilaterals of diameter h K with outward unit normals n#- The element 
boundaries dK have an inflow boundary dK- = {x G dK : (3 ■ < 0} and an outflow boundary 

dK + = dK \ dK The space of admissible solutions is extended to the partition using the broken 
space V(Vh) = HK^p h V(K). The standard conventions in finite element meshing are assumed to 
be in force: Vh is a family of partitions Vh and each element K of Vh is the image of an invertible 
map F k of a master element K = [— 1, l] 2 . The partitions Vh £ Vh are regular and, in the present 
study, it is sufficient to take F K as affine maps. For each partition Vh-, approximate solutions are 
sought in the subspace V p (Vh ) = {v € L 2 (Cl) \ v\ K o F~* G Q Pk (K)} where Q* K (K) denotes the 
space of functions formed by tensor products of Legendre polynomials of degree p K on the master 
element K. Note that the polynomial degree, p K , may vary over different elements in the mesh 
and that functions v p h G V p (Vh) are discontinuous across element interfaces. The approximation 
properties of such spaces are typified by local interpolation estimates of the following type (see [5]): 
if u G H $ (K), there exists a constant C, independent of h K = diam (K) and p K (the minimal order 
of the polynomial shape functions for K), and a polynomial w of degree p K , such that 


i,min(p jr +l,s) 

\\u-w\\ r , K <C - — ; r = 0,1 


K 


( 4 ) 


where ]| • denotes the usual Sobolev norm. 


The following notation is used for functions v G V(Vh)'- 


v 


v ± = 
.int k _ 


,ext k _ 


v 

(v,w) K 

{v,w) dK 


lime -,0 v(x ± e/3) 

u| R .(x), x G dK 

v\ L (x), x G dK H dL 

4- vwdx ; ||u.|| K = y(u,u) K 

IdK vw\(3- n K \ds ; ({v)) dK = \JJy^v)a 


K 


( 5 ) 


The discontinuous Galerkin method applied to (2) is written in the following abstract form: 


Find u G Vp(Vh ) such that 


E B k (u,v) = E . for every t> £ V r (V h ) (6) 

Kev h Kev h 
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where (see [4]) 

B k {u, v) = f (u /3 + au,v+ 8-j-vp) K + (1 + 6^f-){& + - u~, v + ) &K -\t- 


+ (1 + 6-^)(« + ,u + )aif_nr_ 

Pk 

(7) 

- f (/> $ + h-^-vg) K + (1 + 6-f -){g, v)dK-nr- 
Pk Pk 

(8) 


and <5 is a parameter with a value of 0 or 1. The method with 6 = 1 in (7) and (8) is the so- 
called streamline-upwind discontinuous Galerkin method. The additional term tyf-vg in the element 
integrals adds diffusion in the streamline direction without compromising the accuracy of the ap- 
proximation. The method with 6 = 0 is the standard discontinuous Galerkin method which can be 
viewed as a higher-order extension of a cell-centered finite volume method where the coefficients of 
the higher-order terms in the polynomial approximation of the solution in an element are obtained 
from the conservation law and not by reconstruction. Integrating the first-order terms by parts in 
(6) with 6 = 0 and manipulating the result yields the familiar numerical flux formulation of the 
finite volume method 


where 




B k (u, v) — (u, av — vp) K + i ext K )i>ds 

, w ext K ) = \ (i int K + » ext K ) - \\P -n K \ (4 ext K - fi int K 


(9) 

( 10 ) 


The error in the discontinuous Galerkin solution satisfies the following a priori estimate [4]: 


Theorem 1 Let u £ H S {?1) be a solution to (2), let u be a solution to (6), and let (4) hold. Then 
there exists a positive constant C, independent ofh K , p K , and u, such that the approximation error, 
e = u — u, satisfies the following estimate 


ll e ll|ftp,/3 ^ c 



( 11 ) 


where p K — min^^ -f 1, s) — 4, v K = s — 1, and 


Wllw = \ E 

K&V h 


hit < 


, 2 IMIi + ||e||^. + ({e + e ))g K _\ T _ + (( e ))dKnda 
L Pk 


( 12 ) 


The a priori estimate (11) establishes convergence of the method and is useful for predicting 
how the error in numerical solutions behaves with //-refinement or p-enrichment, Unfortunately, 
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its usefulness in assessing the accuracy of a given numerical solution is limited since the estimate 
involves unknown constants and the exact solution. 


A POSTERIORI ERROR ESTIMATION 


A posteriori error estimates used here are based on extensions of an element residual method of 
Ainsworth and Oden [6]. Element error indicators are computed by solving a suitably-constructed 
local problem with the element residual as data. These local indicators are used in the adaptive 
strategy to assess the accuracy of the solution in an element. Moreover, they contribute to a 
global error estimate which is accurate enough to provide a reliable assessment of the quality of the 
approximate solution. Detailed derivation of the a posteriori estimate can be found in [4], 


The local problem is constructed to result in an upper bound on the error. Let ip K be the 
solution to the following local problem, 


where 


Bk^KI V k) ~ ^K^Vk) <{^K^ V k) V U € 


K&k^k) - Wk-,/3’ Vv K ) K + a(ij; K ,v K ) K 

Pk 


(13) 


(14) 


and a > 0 is a constant. Note that the local problem differs from the conservation law, in particular, 
it is symmetric and induces a norm on the space V(K). The solution to the local problem, measured 
in the norm, 

( 15 ) 


serves as an element error indicator in the adaptive strategy. The global error estimate is a sum of 

element contributions given by 

^ ( 16 ) 


M\ aU = JY.Uk\\\ u 


K&Vk 


A (K) 


The solution to the local problem (13) provides an upper bound on the global error in the following 
sense [4]: 


Lemma 1 Let if G V(Vh ) be the solution to the following problem: 

E = E B{e,v) Vu G V(V h ) 

Kev h Ken 

Then there exists a positive constant k such that 

11^11^17 — ^lll e llkpi/3 


(17) 


(18) 
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An approximate solution to the local problem (13) in the corresponding norm serves as a local 
error indicator for the element. Since the discontinuous Galerkin solution satisfies the orthogonality 
condition, 

B K (e, v) = 0 VveQ P *(K) (19) 

the error indicator must be approximated with a polynomial of degree p K + o K where a K > 1 in 
order for the discrete local problem to have a non-trivial solution. If a complete polynomial of 
degree p K + a K (on the master element) is used to approximate the solution to the local problem, 
then the discrete local problem requires the solution of a system of order (p K +cr K + l) 2 . This system 
can be fairly large compared to the system of (p K + l) 2 equations used to obtain the approximate 
solution for which we are estimating the error. Since (p K + l) 2 terms on the right hand side of the 
discrete local problem (corresponding to (19)) are zero, a simplification is made by approximating 
the solution to the local problem in the space Q p k +<t k^K) \ Q p x(K). In other words, the solution 
to the local problem is approximated with incomplete polynomials of degree p K + a K by neglecting 
the terms associated with polynomials of degree p K . This simplification results in a system of 
<J K (a K + 2 p K + 2) equations for each element. 

THE hp - ADAPTIVE STRATEGY 

The hp-adaptive strategy used here is an extension of the 3-step strategy developed by Oden, 
Patra, and Feng [7] for a large class of elliptic problems and, in several applications, was shown to 
yield exponential rates of convergence with respect to both CPU time and the number of unknowns. 

The goal of the adaptive strategy is to deliver a solution with the specified level of error in three 
adaptive steps: (1) estimate the error in the solution obtained on an initial mesh (2) construct a 
new mesh using h-refinement of the initial mesh, solve the problem on the new mesh, and estimate 
the error, and (3) enrich the approximation in regions where the solution is smooth by increasing 
the spectral order of the elements in the mesh from step (2), and if necessary, perform h-refinement 
in regions where the solution is of low regularity. If the level of error after step (3) exceeds the 
specified level, it is necessary to repeat steps (2) and (3) until the desired error is attained. 

The /ip- adaptive strategy is based on the assumption that the a posteriori estimate is a rea- 
sonable approximation to the actual error in a particular solution. The a priori estimate (11) and 
some additional assumptions (see [4]) lead to expressions for estimating the local regularity of the 
solution and for predicting the mesh required to reduce the error to the specified level. The entire 
procedure is outlined below. Detailed development of the hp - adaptive strategy can be found in [4]. 

(i) Specify a target normalized error, rj T . The target error is normalized by the solution 
in the same norm. Specify the parameter a to determine the intermediate target error, 
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rjj = arj T . Specify the parameters a s and a D which establish reduction factors for the 
error in smooth and non-smooth regions as described below. Specify the parameters 
p K and v K in the a priori estimate (11). Formally these parameters depend on the 
global regularity of the solution. While there is little theoretical justification, local 
values can be used by computing the rate of convergence of the local error for a uniform 
/i-refinement and p-enrichment of a coarse mesh. 

(ii) Construct an initial mesh Vo containing N(V o) elements. The elements in Vo have 
uniform p K = po and essentially uniform h K « h 0 . Find the approximate solution 
u 0 € V Po (Vo). Estimate the error Oq where 

«o = J t = ,/E IIVvII’ (20) 

y Kev<> y KeTo 

and if; K is the solution to the local problem (13). 

(iii) Construct a mesh Vi by subdividing each element in Vo into the number of elements, 
n K i required to equally distribute the error and reduce it to 9i = f?j(||«o|l A[r + #o)- The 
number of elements, n K , is obtained by iteratively solving the following two equations: 



N(Vi) = E "k (22) 

Kev o 


Find the approximate solution u\ G V Po (Vi) and estimate the error 0\. 

(iv) Estimate the local regularity of the solution by computing the rate of convergence of 
the local error 


Vk 


logflo,A--log 


log h K - log 


K = l,---,N(Po) 




(23) 


The value of p K given by (23) is associated with an element K in the initial mesh and 
is simply inherited by the new elements generated by subdividing the element K. The 
expected rate of convergence for smooth solutions is p K + according to the a priori 
estimate (11). Divide the error into two contributions according to the value of p K : 


<? D 

Os 


I o * \k ; 

K£U d 

■> 


= {K € Vi : p K < Pk + 

n s = Vi\n D 


(24) 

(25) 
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Subdivide the elements in 0 D into the number of elements required to equally distribute 
the error and reduce it to a D 0 D . Enrich the approximation in Q s to equally distribute 
the error and reduce it to a s 9 s by increasing p K according to 


Pk =Po 


\ a s e s ) 


(26) 


Find the approximate solution on the new mesh and estimate the error. 


(v) If the estimated error in (iv) is larger than the target error, repeat step (iii) and 
until the target error is reached. 


(iv) 


In the current implementation, h-refinement is accomplished by successive bisection of an ele- 
ment and is limited to two levels for a particular adaptive step. The h-refinement in (iv) is necessary 
only when the error 0 D exceeds the target error. 


NUMERICAL EXAMPLES 


The discontinuous Galerkin method with <5 = 1 in (6) is used to solve the model problem (2) 
to assess the reliability of the error estimate and to investigate the performance of the /ip- adaptive 
strategy for problems with discontinuous solutions. 


Example 1 


We solve the linear model problem (2) with the following data: 


(i) 

n = 

(-1,1) X (-1, 

1) 

(ii) 

P = 

(1.0,0.0) T 


(iii) 

a(x) 

= 1.0 




/ 3 g— 5(1+J/ 2 ) 

if y < 0 

(iv) 

g = 

} _3 e -5(l+j/ 2 ) 

otherwise 


The source term / in (2) 


is chosen so that the exact solution is the discontinuous function given by 


u(x,y) 


3e -5(* 2 +y 2 ) if y <0 
— 3e _5 ( a;2+ ^ 2 ^ otherwise 


(27) 
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and shown in Fig. 1. The discontinuity is aligned with element interfaces at y = 0 to illustrate the 
advantage of using a discontinuous method to capture discontinuities, particularly if the adaptive 
scheme includes some shock fitting which aligns the grid with the discontinuity. 

The problem was solved using a variety of uniform meshes with h-refinements, p-enrichments, 
and the /ip- adaptive strategy. The error history for an /ip-adaptive solution starting from an initial 
8x8 mesh of p = 1 elements is listed in Table 1. The target error was nearly achieved at each step 
in the adaptive process. Recall that the target error is a global quantity obtained as the square 
root of the sum of the squares of the element error indicators (16). Therefore the error in a single 
element cannot exceed the target error. The global effectivity index, the ratio of the estimated error 
to the actual error, is also listed in Table 1. An effectivity index close to unity indicates that the 
error estimate is reliable and provides a good approximation to the actual error. The effectivity 
indices in Table 1 are slightly less than unity, indicating that the actual error is larger than the 
estimated error. However, the estimated error is sufficiently close to the actual error to result in an 
effective adaptive strategy. 


Adaptive step 

Target error 

Achieved error 

Effectivity index 

initial (8x8 mesh, p = 1) 

— 

15.4% 

0.998 

h-refinement 

7.5% 

3.3% 

0.996 

p-enrichment 

5.0% 

5.5% 

0.901 


Table 1: Example 1 - Error history for an hp - adaptive solution 


The rate of convergence of the estimated and exact error is compared in Fig. 2. The exact error 
(denoted by a solid line in the figure) and the estimated error (denoted by a dashed line) are in 
close agreement, indicating the reliability of the estimate. Note that with the discontinuity aligned 
with element interfaces, the error behaves as if the solution is smooth; that is, algebraic rates of 
convergence are achieved with respect to mesh refinement, and exponential rates of convergence 
are achieved with respect to p-enrichment. When the discontinuity is aligned with the element 
interfaces, the most significant error reduction with fewest degrees of freedom results by specifying 
a target error for the h-step which is closer to the initial error than to the final target error. This is 
verified by the curves corresponding to two hp - adaptive solutions in Fig. 2. The error corresponding 
to the /ip- adaptive solutions in Fig. 2 exhibits super-linear rates of convergence. 

The element residual method gives a global error estimate which bounds the actual global error, 
however, nothing in the theory indicates the reliability of the local element indicators. Since the 
local error indicators are used in the adaptive strategy, it is important that the indicators give an 
accurate approximation of the actual element error. The hp - adapted mesh resulting from an initial 
8x8 mesh of p = 1 elements and the local effectivity index for the element error indicators, rj K , are 
shown in Fig. 3. Although there are some elements with low effectivity indices (indicating that that 
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actual error is much larger than the estimate), the local effectivity index for most of the elements 
falls between 0.8 and 1.2 indicating that the local error indicators are sufficiently accurate for use 
in the adaptive strategy. 


Example 2 

The following data is used in (2): 


(i) 

ft = 

(- 1 , 1 ) x (' 

-i,i) 

(ii) 




(iii) 

a(x) 

= 1.0 


(iv) 

9(x, 

f 5e” 

li+y 2 ] . 

- { -i 

— 8e~' 


. 1-121 


X — — 1 

y = -l 


The source term / in (2) is chosen so that the exact solution is a function which is discontinuous 
along the domain diagonal given by 


f 5e 1^+2 ) 2 +^ 2 ] -|- 3e J) 2 ] if y>x 
{ — 1 — Se~ 5 ^ x ~^ 2 +( ~ y+ 2 ^ otherwise 


(28) 


and shown in Fig. 4. 

The global estimated error for a sequence of uniform refinements and for several adaptive hp- 
meshes is shown in Fig. 5. The labels hp- adaptive in Fig. 5 refer to the adaptive strategy with 
only p-enrichment in the third adaptive step. The labels hhp - adaptive refer to the strategy with 
both h-refinement and p-enrichment in the third adaptive step. The hp - adaptive strategy delivers 
nearly linear rates of convergence with respect to the number of degrees of freedom. The rates of 
convergence (the slope of the lines in Fig. 5) for the adaptive strategy are higher than the rates of 
convergence for uniform refinement, indicating that a more accurate solution is obtained with far 
fewer degrees of freedom when using the /ip-strategy. The rate of convergence obtained with the 
adaptive strategy determines the efficiency of the overall process, and as seen in Fig. 5, the rate of 
convergence depends significantly on the target and intermediate error specified. 

The error history for the /tp- adaptive solution denoted by the solid triangles in Fig. 5 is listed in 
Table 2. The target error was nearly achieved at each step in the adaptive process. The effectivity 
index (the ratio of the estimated error to the exact error) is on the order of 0.6, quite good for a 
discontinuous solution, but indicating that the actual error is larger than the estimated error. • 
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Adaptive step 

Target error 

Achieved error 

Effectivity index 

initial (16 x 16 mesh, p = 1) 

— 

7.22% 

0.62 

h-step 

3.6% 

4.1% 

0.58 

hp- step 

2.2% 

2.8% 

0.57 


Table 2: Example 2 - Error history for an hp- adaptive solution 


Recall that the global error is a sum of element error indicators. The primary source of the 
under-estimation of the global error is the under-estimation of the element error indicators near 
the discontinuity as shown in Fig. 6. Although the local error estimate provides a qualitative 
measure of the error at the discontinuity, the low local effectivity index indicates some severe 
under-estimation of the error in that region. Note, however, that the local error estimate in smooth 
regions is very accurate with effectivity indices near unity. 


CONCLUDING REMARKS 


The development of an hp- adaptive discontinuous Galerkin method for hyperbolic conservation 
laws is presented in this work. The emphasis of the work is on a model class of linear hyperbolic 
conservation laws for which it is possible to develop a priori error estimates and reliable a posteriori 
estimates which provide bounds on the actual error. These estimates are obtained using a mesh- 
dependent norm which reflects the dependence of the error on the local element size and the local 
order of the approximation. 

The hp- adaptive strategy is designed to deliver solutions to a specified error level in an efficient 
way. This is accomplished using a three-step procedure in which the a posteriori estimate is used to 
determine the error in the solution at -a particular adaptive step and the a priori estimate is used 
to predict the mesh required to deliver a solution with the specified level of error. The hp- adaptive 
strategy makes further use of the a priori estimate to provide detection of discontinuities in the 
solution thereby identifying regions where h-refinement and p-enrichment are appropriate. 

Numerical experiments demonstrate the effectiveness of the a posteriori estimates in providing 
reliable estimates of the actual error in the numerical solution. Although local error estimates 
near discontinuities under-estimate the actual error, the local error estimates are very accurate in 
smooth regions. The numerical examples also illustrate the ability of the hp- adaptive strategy to 
deliver a final solution with the specified error. While the hp-adaptive strategy provides super-linear 
convergence rates with respect to the number of unknowns in the problem, the rate of convergence 
depends on the level of error requested at each step in the adaptive process. More numerical 
experiments are needed to provide guidelines for selecting the optimum user-specified parameters. 
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Figure 5: Example 2 - Rates of convergence of the estimated global error with the number of 
unknowns. 



Figure 6: Example 2 - Reliability of the local error indicators for an hp-adapted mesh. 
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